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Abstract. - We use a lattice Boltzmann (LB) kinetic scheme for modelling amphiphilic fluids 
that correctly predicts rheological effects in flow. No macroscopic parameters are included 
in the model. Instead, three-dimensional hydrodynamic and rheological effects are emergent 
from the underlying particulate conservation laws and interactions. We report evidence of 
shear thinning and viscoelastic flow for a self-assembled gyroid mesophase. This purely kinetic 
approach is of general importance for the modelling and simulation of complex fluid flows in 
situations when rheological properties cannot be predicted a priori. 



Introduction. - Lattice Boltzmann (LB) modelling schemes have emerged in the last 
several years as a powerful approach for simulating the dynamics of a variety of complex 
systems, from flow with suspended particles to multiphase flow through porous media [1]. Of 
great interest is the extension of this method to the modelling of flow properties of viscoelastic 
materials such as polymer melts and amphiphilic fluids, where the dynamics of microscopic 
structures is coupled to the flow. 




Fig. 1 - a) High-density volume rendering of one of the two immiscible fluids (red) and interface 
between the two (blue) in the gyroid mesophase. b) Same as in a), but after 8000 simulation timesteps 
with an imposed steady shear {U = 0.1, a; = in eq. JSJ)- The shear reduces the crystallinity and the 
material assumes more fluid-like properties. Model size is 64^ 
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Qian and Deng [2] correctly described transverse wave propagation using a lattice BGK 
model with an ad hoc modified equilibrium distribution while Ispolatov and Grant [3] ob- 
tained viscoelastic effects by adding a force to represent memory effects in the LB equation. 
However, both approaches involve the inclusion of macroscopic parameters such as Young's 
modulus [2] or elastic coefficients [3]. These methods cannot be regarded as fully mesoscopic 
as at least one parameter is imposed on the basis of macroscopic considerations. Giraud et 
al. [4] introduced a single fluid, two dimensional model to treat a macroscopically viscoelastic 
fluid and later Lallcmand et al. [5] extended this model to three dimensions, but to the best 
of our knowledge no generalization for multiphase flow or flow with suspended particles has 
yet been implemented. A free energy Ginzburg-Landau (GL) model can be defined to study 
rheological properties of complex fluids [6,7]. A popular LB scheme based on the same GL ap- 
proach proceeds by deflning an equilibrium distribution through the imposition of constraints 
on macroscopic thermomechanical quantities such as the stress tensor. Using such a scheme, 
Denninston et al. [8] obtained non-Newtonian flow behaviour (including shear thinning and 
banding) using an LB algorithm for the hydrodynamics of liquid crystals in the isotropic and 
nematic phases. However, with such methods the dynamics is not dictated by the mesoscopic 
processes; numerical instability [9] can make this approach unsuitable for a fully mesoscopic 
description of the dynamics. 

Gonzalez and Coveney [10], using a fully mesoscopic, kinetic approach which does not re- 
quire the existence of a thermodynamic potential obtained a self assembled gyroid mesophase 
in the course of simulating an amphiphilic fluid formed by two immiscible fluids and a surfac- 
tant species (fig.0. The term amphiphilic fluid is used to describe a fluid in which at least one 
species is made of surfactant molecules. Surfactants are molecules comprised of a hydrophilic 
(water-loving) head group and a hydrophobic (oil-loving) tail. An amphiphilic fluid may con- 
tain oil, water, or both fluids in addition to surfactant. Such complex fluids can self-assemble 
to form ordered structures such as lamellae, micelles, sponge and liquid crystalline (cubic) 
phases showing pronounced rheological properties [11]. 

In this letter we use this purely kinetic LB method to model complex flows whose rheo- 
logical properties are emergent from the mesoscopic kinetic processes without any imposed 
macroscopic constraints [12]. In particular, we show evidence of the appearance of rheological 
effects, such as shear thinning and viscoelasticity, for a self-assembled gyroid liquid crystalline 
cubic mesophase. 

The model. - A standard LB system involving multiple species is usually represented by 
a set of equations [13] 

<(x + c„t+l)-<(x,i) = -^«(x,<)-nr«(x,i)), i = 0,l,...,6, (1) 

where nf(x,t) is the single-particle distribution function, indicating the density of species a 
(for example, oil, water or amphiphile), having velocity c,;, at site x on a D-dimensional lattice 
of coordination number b, at time-step t. The collision operator Of represents the change in 
the single-particle distribution function due to the collisions. We choose a single relaxation 
time Tq, 'BGK' form [14] for the collision operator In the limit of low Mach numbers, the 
LB equations correspond to a solution of the Navier-Stokes equation for isothermal, quasi- 
incompressible fluid flow whose implementation can efficiently exploit parallel computers, as 
the dynamics at a point requires only information about quantities at nearest neighbour lattice 
sites. The local equilibrium distribution n""^"^ plays a fundamental role in the dynamics of the 
system as shown by eq. In this study, we use a purely kinetic approach, for which the 
local equilibrium distribution n"'^'' (x, t) is derived by imposing certain restrictions on the 
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microscopic processes, such as explicit mass and total momentum conservation [15] 



2ct 



u 
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(2) 



where u = u(x, i) is the macroscopic bulk velocity of the fluid, defined as n"(x, i)u" = 

(x, i)ci, are the coefficients resulting from the velocity space discretization and 
is the speed of sound, both of which are determined by the choice of the lattice, which is 
D3Q19 in our implementation. Immiscibility of species a is introduced in the model following 
Shan and Chen [16,17]. Only nearest neighbour interactions among the immiscible species are 
considered. These interactions are modelled as a self-consistently generated mean field body 
force 



F"(x,t) = -^"(x,t) 



.^V^(x',i)(x'-x) 



(3) 



where ip"{yi,t) is the so-called effective mass, which can have a general form for modelling 
various types of fluids (we use tp" = {1 — e~" ) [16]), and gsa is a force couphng constant 
whose magnitude controls the strength of the interaction between components a, a and is set 
positive to mimic repulsion. The dynamical effect of the force is realized in the BGK collision 
operator in eq. (??) by adding to the velocity u in the equilibrium distribution of cq. |(5Jl an 
increment 



As described above, an amphiphile usually possesses two different fragments, each having 
an affinity for one of the two immiscible components. The addition of an amphiphile is 
implemented as in [12]. An average dipole vector d(x, t) is introduced at each site x to 
represent the orientation of any amphiphile present there. The direction of this dipole vector 
is allowed to vary continuously and no information is specified for each velocity , for reasons 
of computational efficiency and simplicity. Full details of the model can be found in [12] 
and [18]. 

Simulation details. - We use LB3D [19], a highly scalable parallel LB code, to implement 
the model. A single simulation of a 64"^ model, i.e. a point in fig[5] below, needs around 300 
Mbytes of RAM and takes about 50 CPU hours to complete on a single processor machine. The 
very good scaling of our LB3D code permits us to run all our simulations on multiprocessor 
machines and computational grids in order to reduce the length of data collection to a few 
weeks. The simulation parameters are those used in [10] that lead to the formation of gyroid 
mesophases. We use a 64^ lattice size with periodic boundary conditions. The simulations 
were all started from a checkpointed configuration of a 200000 timesteps equilibrated gyroid 
mesophase [19]. A single unit cell of the gyroid minimal surface is of the order of 5 — 8 lattice 
lengths in its linear dimensions. This implies that hundreds of unit cells are present in a 64^ 
sample. Therefore, finite size effects are not expected to affect the qualitative outcome of 
the simulations. Moreover, some simulations with a 128'^ lattice have been run to confirm 
the results obtained and the 128"^ simulations reported in [20] are consistent with the results 
shown here. Nevertheless, such lattice sizes oblige us to study perfect gyroids; much larger 
system sizes (and computational cost) would be necessary in order to investigate the role of 
domains and their associated defects [21]. We leave this for future work. 

In order to inspect the rheological behaviour of gyroid mesophases, we have implemented 
Lees-Edwards boundary conditions, which reduce finite size effects if compared to moving solid 
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Fig. 2 - Plot of Oxz (dimensionless units) vs. shear rate ^xzTr (Weissenberg number) for an LB 
simulation of a 64^ gyroid mesophase. is the stress relaxation time; see text for discussion. The 
inset shows the percentage drop of the viscosity as the shear rate "fxzT,- increases (shear-thinning). 

Fig. 3 - Plot of axz (black dots) vs. time for oscillatory shear [U = 0.05 and u = 0.0004, 0.005 in 
eq. jnj) applied to 64'^ lattice-Boltzmann gyroid mesophase. All quantities are dimensionless. After 
a brief transient, Oxz oscillates with the imposed frequency ujTr- The strain (dashed line) is plotted 
(scaled) as a guide to the eye. The continuous line is the best fit line for the a^z data. Shear moduli 
G'{u!rr), G" (ujTr) are calculated as in eq. Q over a time interval shown by the continuous horizontal 
line. 

walls [22]. This computationally convenient method imposes new positions and velocities on 
particles leaving the simulation box in the direction perpendicular to the imposed shear strain 
while leaving the other coordinates unchanged. Choosing z as the direction of shear and x as 
the direction of the velocity gradient, we have 

{{z + A^) mod ,x> ( u,^ + U ,x> 

z mod N;, ,0 < X < u',_ = I ,0<x<N^ , (5) 

(z-A^)mod7V3 ,x<0 [ u,_ -U ,x <0 

where = UAt, U is the shear velocity, is the z— component of u and iV^-j-j.) is the system 
length in the x{z) direction. We also use an interpolation scheme suggested by Wagner and 
Pagonabarraga [23] as A^ is not generally a multiple of the lattice site. Gates et al. [24], found 
pronounced artefacts (lock-ins) in simulations of 2D sheared binary mixtures with multiple 
Lees-Edwards planes. In our own work, we have not seen any of these artefacts, even in the 
longest simulations performed, and a linear velocity profile is obtained at steady state. 

Consistent with the hypothesis of the LB model, we set the maximum shear velocity to 
U = 0.1 lattice units. This results in a maximum shear rate jxz = ^ g"' = 3.2 x 10^"^ in 
lattice units. Simulations are run for T — 10000 timcsteps in the case of steady shear. Steady 
state is reached in the first 1000 timcsteps, and the relevant component of the stress tensor 
is averaged over the last 3000 timcsteps. Some measurements were repeated by doubling the 
simulation time to a total of T = 20000 timcsteps but no significant differences were found. 
For oscillatory shear, we set 

U{t) = Ucos{Lut) , (6) 

where ti;/27r is the frequency of oscillation and U — 0.05 lattice units. These simulations were 
run for at least three complete oscillations, t ~ 5000 — 100000 timcsteps. 

Results. - In order to investigate rheological properties of the system, we perform simu- 
lations of flow under shear, mimicking a rheometer. We set U = n x 0.005, n = 2 ... 21 and 
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a; = to impose a stationary Couette flow. Once steady state flow is reached, we collect the 
relevant component of the pressure tensor Uxz- According to Newton's law for viscous flow of 
a liquid we have 

Cxz = 2?77a;2 , (7) 

where is the viscosity of the liquid and = i^slisi^iliil . In fig.|21we plot ct^z against the 
imposed non-dimensionalized (see text below) shear rate '^xzTr (Weissenberg number). We 
note that the slope of the curve changes as these shear rate increases. This indicates that the 
viscosity 77 depends on the shear rate '^xz, V = viixz)- We therefore conclude that our model 
exhibits non-Newtonian flow behaviour. Complex fluids such as amphiphilic mixtures are well 
known to exhibit such behaviour experimentally [11]. In the inset we calculate the percentage 
drop of 77 by averaging the slope of the curve over subsets of four data points. We note that 
the viscosity decreases as the shear rate increases. This behaviour is referred to as "shear 
thinning" [11]. In general, fluids exhibit macroscopic non-Newtonian properties because of 
underlying complex mesoscopic interactions and leading to changes in their microstructural 
properties. In our case, the highly ordered structure of the gyroid mesophase is responsible 
for this rheological effect, and our LB model correctly captures it. Therefore, considering the 
underlying physics of mesophases, this simulation result is not only justified but also expected. 
We note that the predictions of the model could be directly verified, but no experimental 
evidence for gyroid rheology is currently available in the literature. 

For oscillatory shear, in eq. © we set U = 0.05 lattice units and we span a range of two 
decades of frequencies by varying cu between u; = 0.0001 and uj — 0.01 lattice units. For a 
viscoelastic medium 

a^z = 7o[G'(w) sm{Lut) + G"{uj) cos(tjt)] , (8) 

where 70 = is the imposed shear strain and G'{uj), G"{uj) are, respectively, the storage 
and loss moduli which respectively measure the elastic and viscous response at any given 
frequency [11]. In fig. we show the time dependence of a^z for two simulations with different 
values of uj. After a brief transient, axz begins to oscillate as predicted by eq. l(Hl with the 
imposed frequency u and an amplitude ctq and shift (j) that depend on the storage and loss 
moduli 



^, ^ gp cos(0) ^„ ^ tjQ sinjcj)) 
7 7 ' 

where uo and <j) are the fitted values for the amplitude and phase shift. We derive G'{uj) 
and G"{uj) for the two different values oi uj ~ 0.005,0.0004 by fitting ctq and (p with the 
standard Levenberg-Marquardt algorithm over at least three decades (see caption for details) . 
In fig. ^ we plot a^z against the shear strain at different times (Lissajous plots) for the 
two frequencies of figOl We note that for the higher frequency (left panel), axz is in phase 
with the strain, indicating that the gyroid mesophase exhibits a solid-like response at short 
timescales. For the smaller frequency, the phase shift is almost 7r/2, typical of a liquid-like 
response as dictated by eq. Q. In fig. |31 we plot G'{ujTr) and G"{uJTr) for the different 
simulation values of uj. In order to capture the relevant scales in the model, all quantities are 
dimensionless. In particular, the frequency uj has been multiplied by the characteristic stress 
relaxation time, ~ 220 simulation timesteps, calculated by fitting the exponential decay of 
axz when shear is removed. We note that a crossover between G'{ujTr) and G"{uJTr) occurs 
at ujTr ~ 0.05. This signals a transition between solid and liquid-like behaviour and indicates 
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Fig. 4 Fig. 5 

Fig. 4 - Lissajous plots of a^z (in lattice units) vs. dimensionless strain ^^z- In the left(riglit) panel, 
points tend to form a straight line (ellipse, respectively). This indicates that the pressure a^z is in(out 
of) phase with respect to the strain 'y^z, as typical for solid(liquid)-like phases. The overall behaviour 
is typical of viscoelastic fluids. 

Fig. 5 - Plots of shear moduli G' (ujTr) (regular triangles), G" {ujTr) (inverted triangles) calculated 
for all frequencies used in our simulations. We note that a crossover between liquid and solid-like 
behaviour occurs around ui = 0.5. Theoretically, G' ~ a;'^, G" ~ w as goes to 0. Our data agree 
with this prediction as shown by the dashed lines in the plot. All quantities are dimensionless (Go is 
the plateau modulus for G'). 



that our LB model is capable of predicting viscoelastic flow; see also fig. Q and caption. We 
note that a dimensionless crossover point of wr^ ~ 1 consistently links rheological properties 
of the mixture with the relaxation time of the gyroid mesophase. This is encouraging for the 
use of this simulation model to inspect the correlation between macroscopic dynamics and 
mesoscopic structure. Indeed in the case of linear viscoelasticity, G' ~ w^, G" ^ u) as lo goes 
to for any fluids. In fig. |Slwe plot lines ~ lo and ~ as a guide for the eye; our data show 
excellent quantitative agreement with theory. We interpret the elastic response component as 
due to the stress response to mechanical perturbations of the long-range ordered equilibrium 
structure of the liquid gyroid crystalline mesophase. We also tried to fit the data in figEl to 
a simple single relaxation time (Maxwell) model of viscoelasticity but the fit is very poor. 
This suggests that there are several relaxation times present, as is to be expected here since 
the viscoelasticity arises from complex mesoscale structure. We note that the existence of a 
spectrum of relaxation times is consistent with the stretched-exponential behaviour of domain 
self-assembly in amphiphilic fluid systems found in [10]. 

Conclusions. - In this letter we use a purely kinetic LB model that leads to the emergence 
of non-trivial rheological properties. Non-Newtonian behaviour, in this case a decrease of 
viscosity with shear rate, has been shown for an initially self-assembled gyroid liquid crystalline 
mesophase. In addition, linear viscoelastic effects in the system are manifest in the simulations. 
It is notable that our model correctly predicts the theoretical limits for the moduli G'{lo), 
G"{llj) as uj goes to as well as a crossover in G'{uj), G"{uo) at higher lo. We note that, unlike 
previous approaches, this model does not require any assumptions at the macroscopic level, 
that is it provides a purely kinetic-theoretical approach to the description of complex fluids. 
This model can therefore help in understanding complex flows, such as flow of viscoelastic 
liquids in porous media, colloidal fluids and polymer melts. Work is in progress to assess 
the generality of this approach within our lattice-Boltzmann model of amphiphilic fluids. 
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In particular, such simulations should enable us to investigate the link between mesoscopic 
structure and macroscopic dynamics, for example by correlating rheological relaxation times 
to mesoscale relaxation processes within the amphiphilic fluid. 
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